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We present a detailed study of the first simple mechanical system that shows fully realistic trans- 
port behavior while still being exactly solvable at the level of equilibrium statistical mechanics. The 
system under consideration is a Lorentz gas with fixed freely-rotating circular scatterers interacting 
with point particles via perfectly rough collisions. Upon imposing a temperature and/or a chemical 
potential gradient, a stationary state is attained for which local thermal equilibrium holds for low 
values of the imposed gradients. Transport in this system is normal, in the sense that the trans- 
port coefficients which characterize the flow of heat and matter are finite in the thermodynamic 
limit. Moreover, the two flows are non-trivially coupled, satisfying Onsager's reciprocity relations to 
within numerical accuracy as well as the Green-Kubo relations . We further show numerically that 
an applied electric field causes the same currents as the corresponding chemical potential gradient 
in first order of the applied field. Puzzling discrepancies in higher order effects (Joule heating) are 
also observed. Finally, the role of entropy production in this purely Hamiltonian system is shortly 
discussed. 

PACS numbers: 05.60.-k, 05.20.-y, 44.10. +i 



I. INTRODUCTION 

One of the major motors of research in statistical 
physics has been the quest to understand the link be- 
tween macroscopic phenomena and the underlying mi- 
croscopic physics of a system. In particular, the origin 
of the macroscopic "laws" of thermodynamic transport 
is still one of the major challenges to theoretical physics. 
These phenomenological laws are known to describe ac- 
curately the processes of diffusion, heat conduction, vis- 
cosity among a host of other phenomena, and are fun- 
damental to the quantitative description of macroscopic 
systems in general. However, the attempts to link these 
macroscopic laws to the underlying microscopic dynam- 
ics have not been conclusive thus far. From a mathe- 
matically rigorous point of view, very few results have 
been obtained Indeed, to our knowledge, the validity 
of Fourier's law has been proven analytically only for a 
very specific model in the limit of infinite dilution with fi- 
nite mean free path ^ . There have also been attempts 
to link transport phenomena to the chaotic properties of 
the underlying classical dynamics [0; and a connection 
between the rate of entropy production and the rate of 
contraction of phase space volume in thermostated (non 
Hamiltonian) systems has been pointed out || ^, 0, g] . 

Given this state of affairs, a common strategy is to 
propose and study systems which reproduce, in numeri- 
cal simulations, the phenomena under consideration, and 
to attempt to determine how the physical ingredients 
of these models give rise to the macroscopic behaviour. 
However, the systems considered thus far have been ei- 
ther too complicated to shed much light upon the prob- 
lem, or have actually failed to reproduce the macroscopic 
phenomenology. Attempts have been made, on one side, 
through the simulation of realistic many-body systems 
satisfying a thermostated dynamics (see for example ||). 



These simulations have indeed been able to reproduce 
non-trivial transport phenomena. However, such studies 
do not provide a detailed understanding of the micro- 
scopic processes involved due to the excessive complexity 
of the system under study. The other numerical approach 
involves the study of transport in "simple systems" . Ex- 
amples of these include: chains of anharmonic oscillators 
pOf and the so called ding-a-ling and ding-dong mod- 
els |ll], ^2|, among others. Of these, energy transport 
in the chains was shown to be "anomalous" . This is 
an euphemism indicating that energy transport cannot 
be described as a diffusive process and currents do not 
scale with the gradients in the expected way. In contrast, 
"normal" transport indicates that even in the thermody- 
namic limit, transport can accurately be described as a 
diffusive process and that the flux is proportional to the 
gradient as stated, say, in Fourier's law. The ding-a-ling 
and ding-dong models, indeed do yield normal transport 
under certain conditions; however, they include geomet- 
ric constraints that make even their equilibrium proper- 
ties an extremely complicated affair. As we shall see, the 
knowledge of such equilibrium properties is often useful, 
which is why we consider it important to have a model 
where these are explicitly known. 

One particularly thorny problem in attempting to re- 
produce the phenomenological laws of thermodynamic 
transport, is that these apply to systems which are 
characterized by local thermodynamic equilibrium (LTE) 
||, Q. That is, systems that are not in equilibrium, 
but for which the intensive thermodynamic variables are 
well defined at each point of the system, and the rela- 
tions amongst these variables are the same as in equilib- 
rium thermodynamics. Thus, results have been obtained 
concerning "diffusive" (normal) energy transport in the 
Lorentz gas |l5| . Yet this system is not described by LTE 
Jl6| . It is therefore not clear what precise meaning can be 
attached to the local temperature appearing in Fourier's 
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law in this situation. 

Under these circumstances, the minimal ingredients re- 
quired in the microscopic physics of a simple system to at- 
tain normal thermodynamic transport in low dimensions 
are presently still under discussion (l], [l^, However, 
we believe that ideal candidate models should possess the 
following features: 

1 . The microscopic physics should be defined in terms 
of reversible Hamiltonian dynamics, since this is the 
nature of known fundamental processes. In partic- 
ular, reversible systems showing an average rate of 
phase space contraction, while of considerable in- 
terest for simulating transport, do not provide a 
fundamental microscopic model. 

2. The equilibrium properties of the model should be 
well understood. 

3. When driven weakly out of equilibrium, they must 
be consistent with the hypotheses of LTE and, of 
course, they must give rise to realistic macroscopic 
transport. 

In this work we present a detailed description and 
extensive simulations of the equilibrium and transport 
properties of a simple reversible Hamiltonian model sys- 
tem which we believe is an ideal candidate to begin to 
unravel the puzzle of thermodynamic transport. The first 
results concerning some of the transport features of this 
model were presented in Jl9[ | . In Section [n| we present 
a detailed description and discussion of the model. We 
show that its equilibrium statistical mechanics is simply 
that of an ideal gas. 

In order to carry out the equilibrium simulations, as 
well as to study the transport properties of the model 
when the system is driven out of equilibrium, it is neces- 
sary to couple the system to thermo-chemical baths. We 
describe the model baths we have implemented to do this 
in Section pll 

In Section|IV| we corroborate through simulations that 
the equilibrium state of the system in the three canonical 
ensembles is indeed an ideal gas. We also present nu- 
merical evidence to show that, when subjected to a weak 
temperature and/or chemical potential gradient, our sys- 
tems reaches a non-equilibrium steady state (NSS) which 
is characterized by LTE. In this situation, our model 
supports both heat and matter flows. We show in Sec- 
tion |y| that these flows are characterized by transport 
coefficients which are independent of system size (that is, 
transport is normal), and that the corresponding trans- 
port coefficients satisfy Onsager's reciprocity relations. 
The fact that the system is homogeneous in energy al- 
lows the prediction of the dependence of the transport 
coefficients on the temperature. We also establish numer- 
ically their dependence on the particle density and dis- 
cuss briefly the effects of subjecting the system to a mag- 
netic field. In Section VI we verify that the Green-Kubo 
formulas connecting transport coefficients with time cor- 
relation functions apply for this system to within numer- 



ical precision, and we compare the effect of an applied 
electric field with that of an applied chemical potential 
gradient, which, to linear order, should induce the same 
flows in the system. In Section VII a brief discussion of 
the applicability of microscopic interpretations of entropy 
production as motor of transport in this system. Finally 
we present a brief summary and mention how this model 
can be modified to study other physical transport prob- 
lems from a microscopic approach. 



II. DEFINITION OF THE MODEL 

In |H| , a Lorentz gas model with elastic collisions was 
used to study heat transport in a quasi-one dimensional 
channel placed between two thermal reservoirs at differ- 
ent nominal temperatures. While the results were consis- 
tent with some sort of "diffusive" energy transport, iden- 
tification with Fourier's law was unfounded: the system 
does not satisfy the hypothesis of LTE and, therefore, 
one cannot define a local temperature, as was shown in 
pq . There, the authors argue that the system does not 
attain LTE due to the existence of an infinite number 
(in the thermodynamic limit) of conserved quantities in 
the dynamics. Indeed, the energy of each particle is con- 
served throughout the evolution of the system. This, in 
turn, implies a breakdown of crgodicity and the result- 
ing process is closer to "colour" diffusion than to heat 
transport. 

The model we study in this work, introduced in |f9| ] , is 
a modification to the usual Lorentz gas model in which 
the scatterers are allowed to exchange energy with a set 
of (non-interacting) point particles of mass m through 
the scattering events. The geometry we consider is that 
of a periodic Lorentz gas in which the hard disc scatterers 
of radius R are fixed on a triangular lattice, the details 
of which are discussed below. The possibility of energy 
exchange is achieved as follows: each disc is a free rota- 
tor with a moment of inertia 0, and scattering proceeds 
according to rules characterizing "perfectly rough" colli- 
sions, that are reversible, conserve total energy and an- 
gular momentum. These collision rules are given by the 
following formulas, which relate the normal and tangen- 
tial components of the particle's velocity v with respect 
to the disc's surface, and the disc's angular velocity uj be- 
fore (unprimed quantities) and after (primed quantities) 
the collision 
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These rules define a deterministic, time-reversible, 
canonical transformation at each collision. The parame- 
ter rj, defined as the ratio between the moment of inertia 
of the disc and the mass of the particle times the square 
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radius of the disc: 
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is the only relevant adimcnsional parameter characteriz- 
ing the collision. It determines the energy transfer be- 
tween discs and particles in a collision. For finite values 
of 77, particles in the system may exchange energy among 
each other through the discs, even though they do not 
interact directly. This simple energy exchange mecha- 
nism overcomes the objections raised in jl(| and permits 
the system to reach thcrmodynamical equilibrium, as we 
will see in Section IV. As r\ — > 0, the rotational energy 



of the scatterers becomes negligible and the collision be- 
comes perfectly elastic, recovering the dynamics of the 
usual Lorentz gas model. For r\ — > 00, the energy of the 
scatterers is unaffected by the collisions with the parti- 
cles, and the angular velocity of the scatterers remains 
constant. In this limit, once again, the system does not 
equilibrate. Thus, in either limit, the energy-mediating 
effect is suppressed and thermodynamical equilibrium is 
not reached. In the following, unless the contrary is ex- 
plicitly stated, we shall always be dealing with the case 
r\ = 1, since that is a value of 77 for which energy ex- 
change, and therefore equilibration, is quite efficient. 

The geometric disposition of the scatterers in the sys- 
tems we study in this work is indicated in Fig. [t| The 
centers of the scatterers arc fixed on a triangular lat- 
tice, along a narrow channel of height 2W, where W is 
the distance between the centers of the scatterers. For 
convenience, this distance was chosen as W = 4R/y/3 
(known in the literature as the critical horizon), which is 
the largest separation for which a particle cannot travel 
arbitrarily large distances without undergoing a collision. 
In this geometry, a set of non-interacting point particles 
of mass m moves freely between collisions with the hard 
discs. In the vertical direction the channel contains two 
discs and periodic boundary conditions are imposed. 

The fact that we put two discs in the vertical direction 
merits comment as, in the simpler "single cell" geometry 
with periodic boundary conditions, spurious effects may 
arise from multiple successive scatterings of one particle 
with the same disc. Some of these effects have been ob- 
served in [ pH , where the same set of collision equations 
has been used to model a deterministic thermostat. We 
have verified that if one considers a system consisting of 
one cell with periodic boundary conditions, containing a 
single scatterer and a single particle, the resulting dy- 
namics gives rise to regular structures in the phase space 
of the particle. In this situation, the system does not 
appear to be ergodic for arbitrary values of 77. We also 
noted that, if instead of periodic boundary conditions 
one considers specular reflections at the boundaries, the 
effect in the particle's trajectory is a "randomization" of 
the tangential component of the particle's velocity in the 
next collision with the disc, recovering a seemingly er- 
godic phase space. Though we believe that the presence 
of other particles will destroy the regular structures that 




FIG. 1: Schematic illustration of the scatterer geometry: the 
scatterers are disposed on a triangular array with finite hori- 
zon to avoid infinitely long trajectories. For matters of con- 
venience, in this work the separation between scatterers has 
been set to have the critical horizon. Periodic boundary con- 
ditions are used in the vertical direction. To avoid spurious 
effects arising from multiple consecutive scatterings of a par- 
ticle off the same disc, we have put two discs on each vertical. 
To study the dependence with system size, the length L of the 
sample is varied. The quoted length is the number of discs. 



appear in the single particle case, we decided to avoid the 
possibility of multiple consecutive collisions of a particle 
with the same disc by placing two discs in the vertical 
direction. 

A central aspect of this model is that its equilibrium 
properties are still trivial, even though, strictly speaking, 
it is an interacting many particle system. Indeed, in spite 
of the modified collision rules, the energy of the system 
is given by 



B = m E >E ^ 
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Thus, all statistical mechanics calculations for this 
system coincide with those of a system consisting of 
a two dimensional ideal gas plus a collection of non- 
interacting free rotors. Hence, a first test for the sys- 
tem is to verify that it equilibrates to a state in which 
its statistical properties are indeed those predicted by 
equilibrium statistical mechanics. In particular, in mi- 
crocanonical simulations, the particle velocities should 
reach Maxwellian distributions at a uniform "temper- 
ature" consistent with the equipartition theorem fl22[ . 
These temperatures should also characterize the distribu- 
tion of angular velocities of the rotating scatterers. The 
same should be true in canonical and grand canonical 
simulations, where the temperature and particle density 
(or, more formally, the chemical potential divided by the 
temperature) are now those established by the values of 
the baths. Furthermore, the equations of state character- 
izing ideal gases should describe the thermodynamics of 
our system. If any of these tests fails, it could be argued 
that the system fails to equilibrate probably due to a lack 
of ergodicity. Unfortunately, passing all the tests does 
not prove that the system is ergodic. We have not suc- 
ceeded in showing rigorously that this system is ergodic. 
However, these numerical tests are fairly stringent: we 
shall see that for the case of an imposed external mag- 
netic field, where ergodicity is known to be violated, an 
effect indeed appears in the energy distribution function 
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of the particles. 



III. THERM O-CHEMICAL BATHS 

In numerical studies of statistical mechanics systems in 
general, and in those on the microscopic origin of macro- 
scopic transport in particular, it is necessary to define 
models for the thcrmodynamical reservoirs in addition to 
the dynamics of the system. The design of these reser- 
voirs and the way these interact with the system, has sev- 
eral subtleties that can produce confusing and, very of- 
ten, wrong results, as has been pointed out in [g3|, gij. In 
this Section, we explain the design of a stochastic model 
that simulates a thermo-chemical bath, which is able to 
exchange energy and particles with the system at fixed 
nominal values for the temperature T and chemical po- 
tential /!. 

The heat and matter reservoir is taken to be an infinite 
ideal gas at temperature T and density p which is placed 
in direct contact with the system being studied. This is 
achieved by assuming that the wall separating the system 
from the reservoir is perfectly transparent to particles 
impinging on it from either side. Of course, it is not nec- 
essary to simulate explicitly the infinite ideal gas which 
acts as the thermodynamical reservoir, instead, we can 
implement this set up using the following rules: when- 
ever a particle of the system impinges on the boundary 
which separates it from the bath, it is removed. On the 
other hand, with a frequency 7, particles are generated 
at the boundary, with a velocity distribution 
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reflecting the assumption that the bath is an ideal gas 
(here and in the following, we take the particle mass and 
Boltzmann's constant equal to one). The choice of T in 
these equations defines the nominal temperature of the 
bath. Moreover, this way of implementing the thermo- 
chemical baths also fixes a nominal value for their chem- 
ical potential. Again, as the bath is an ideal gas, the rate 
7 is given by: 



7 = 



y.1/2 
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Thus, we can express the chemical potential of the bath 
in terms of the parameter 7 as: 
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where Ao is an irrelevant constant. 

Equations (|j), (||) and (||) completely define the algo- 
rithm for the stochastic emission process of the thermo- 
chemical baths used in the simulations, and allow us to 



control the chemical potential and temperature of the 
walls by varying the rate 7 and the temperature. These 
walls were used for grand canonical simulations in equi- 
librium and coupled heat and matter transport simula- 
tions in the NSS. 

To simulate the canonical ensemble in equilibrium and 
pure heat flow in the NSS, diathermal impermeable walls 
are required. These were achieved by reflecting each par- 
ticle that impinged on the wall back into the system, with 
its velocity updated according to the distribution func- 
tion (^) , which again defines the nominal temperature of 
the thermal bath. 

Finally, in order to perform simulations in the micro- 
canonical ensemble, insulating walls are required. We 
have simulated these either by simply considering that 
the particles always perform elastic specular reflection at 
the boundary, or by having no walls at all and considering 
periodic boundary conditions. 

A generalized model for thermo-chemical bath, as well 
as a general approach allowing to verify the validity of 
a given procedure for generating such baths is given in 
Appendix |A|. 



IV. EQUILIBRIUM AND LOCAL THERMAL 
EQUILIBRIUM 

In this section we show that the system described in 
Section [n] reaches a well defined equilibrium state for the 
different equilibrium ensembles. For each ensemble we 
place the system in contact with the appropriate wall, as 
described in the previous Section. However, as expected, 
the equilibrium state that the system reaches does not 
depend on the particular choice of the ensemble. Further- 
more, when the nominal values for the thermodynamical 
quantities fixed by the baths impose a gradient in temper- 
ature and/or in the chemical potential, our model reaches 
a well defined steady state. 

In order to show that our model reaches a satisfactory 
equilibrium state we have measured the velocity distribu- 
tion of the particles P(v) and the angular velocity distri- 
bution of the discs P(u>) in a microcanonical simulation 
in which the particles undergo specular reflections with 
the walls at the left and right extremes of the channel. 
We have fixed the energy of the system and distributed 
it randomly among discs and particles. After some re- 
laxation time the mean energy per particle is twice the 
mean energy per disc, thus satisfying the equipartition 
theorem. In Fig. 0, the measured distributions P(v) 
and P(u>) are shown. The solid line is the correspond- 
ing Boltzmann distributions at the expected temperature 
(T = 150 in arbitrary units). The agreement indicates 
that both discs and particles have reached a state con- 
sistent with the predictions of equilibrium statistical me- 
chanics. This agreement also allows us to relate the av- 
erage energy per particle with the temperature, which 
would have been unfounded otherwise. 

We have computed the particle density and tempera- 
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FIG. 2: Discs angular velocity distribution P[uf) and parti- 
cles velocity distribution P(v) obtained in a microcanonical 
simulation with 30 particles in a channel of length L — 30. 
The solid lines correspond to the Boltzmann distribution at 
the expected temperature T = 150. 
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FIG. 3: Particle's temperature profile T(x) for a canonical 
simulation. The number of particles was N — 30 and L = 30. 
The nominal temperature in both baths was set to same value 
T — 150. The open circles correspond to the time averaged 
temperature of the discs. 



ture profiles in the following way: We divide the channel 
in 1Z disjoint slabs of width Ax = L/1Z. As before, L is 
the length of the channel in the ir-direction. The parti- 
cle density, n(x), in each slab is computed as the time 
average 



n(x)dx 
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x — Xi(t)) dt dx 
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where :Ej(t) is the position of the i-th particle at time 
t and N is the total number of particles in the channel 
that, in the case of a grand canonical situation, depends 
on time. The time average is performed after the steady 
state has been reached. 

Similarly, the time averaged energy density e(x) is 
computed as 
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Here Ei(t) is the energy of the i-th particle at time 
t. From (|?]) we obtain the particle's density profile p(x) 
as the mean number of particles found in the slab which 
contains the position x 



p( x ) — / n(x) dx . 
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The integral in (j^j is taken over the domain of the 
slab that contains position x. Analogously, the average 
particle energy E(x) at the slab containing position x is 
given by 



E{x) = / e(x)dx 

J Ax 
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With (H) and (Q, we calculate the particle tempera- 
ture profile T{x) as 



T(x) = 



p{x) 



(11) 



In Fig. H we show the temperature profile obtained in 
a canonical simulation where both baths were set to the 
same nominal temperature. The temperature obtained 
according to (|ll|) has additionally been averaged over 
an ensemble of 500 different realizations. The number 
of particles was set to N — 30 in a channel of length 
L = 30. We observe that the particles reach an equi- 
librium state characterized by a constant temperature 
along the channel which coincides with the nominal val- 
ues of the baths' temperature. Moreover, the open circles 
in Fig. U correspond to the time averaged energy of the 
discs. The agreement of both profiles indicates equilibra- 
tion between particles and discs. 

In this equilibrium state, the particle's density profile 
(not shown), is also flat. The same behaviour was also 
obtained in grand canonical simulations, as was to be 
expected from the equivalence of the different statistical 
ensembles. 

Let us finally note the following: we have also per- 
formed simulations involving an external magnetic field. 
In this case, it is immediately clear that the system is 
not ergodic. Indeed, there exist isolated circular orbits 
which do not touch any disc. Since particles do not inter- 
act directly, any particle originally on one of these orbits 
will remain for ever. Similarly, this set of orbits cannot 
be reached from initial conditions for which each particle 
touches one disc at least once. Of course, for small fields, 
the orbits that do not touch any scatterer only occur 
when the particle's kinetic energy is low, but since the 
particles can exchange energy with the discs, the kinetic 
energy of a particle can come arbitrarily close to zero, and 
there will be regions that become unreachable for this 
particle. There is thus a true lack of ergodicity for this 
system for arbitrary (non zero) magnetic fields and at all 
temperatures, although the effect becomes weaker as the 
field decreases. We have plotted in Fig. || the distribu- 
tion of particle energies in the microcanonical ensemble 
and clear deviations from the Boltzmann distribution are 
observed. This comes as an indication that the tests we 
have applied should disclose a lack of ergodicity in the 
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FIG. 4: Particle's energy distribution obtained from a micro- 
canonical simulation with an applied external magnetic field. 
The magnitude of the field B is such that qB/mc = 10, where 
q and m are the electric charge and mass of the particles re- 
spectively and c is the speed of light. The lack of low energy 
particles is evident. In the inset, P(-E) is shown on a larger 
energy domain. 



FIG. 5: Temperature profile T(x) with temperature and 
chemical potential gradients. The nominal values for the 
baths' temperature are To = 140 and Ti = 160 with an end- 
to-end chemical potential difference of A(/i/T) = —0.2. The 
mean number of particles was « 25.2 in a channel of length 
L = 30. The open circles correspond to the mean kinetic 
energy of the scatterers. 



system without magnetic field if such were present. 

We now turn to the more interesting out of equilibrium 
situation. In these simulations we connect the two ends 
of the system with two different baths, each of which has 
given values of the chemical potential and the tempera- 
ture. To drive the system out of equilibrium, the nom- 
inal values of the temperatures and chemical potentials 
of these baths are set to differ by fixed amounts. Under 
these conditions, the system is allowed to evolve until a 
non-equilibrium steady state (NSS) develops. In Fig. ||, 
we show results for a typical simulation; the profile cor- 
responds to the average energy per particle for discs and 
particles obtained in a channel of length L = 30. The 
temperature difference of the baths was set to AT = 20 
around a central value T — 150, with a chemical po- 
tential difference of A(/x/T) = -0.2. The profile of the 
average energy per particle is linear and coincides at the 
boundaries with the nominal baths' temperature. As in 
equilibrium, in the steady state discs and particles locally 
equilibrate to the same value of the local mean energy per 
particle, giving a first indication of the establishment of 
LTE. As we will show further on, we will be justified in 
identifying the mean kinetic energy per particle with the 
local temperature. 

We now turn to the verification that LTE holds for our 
system. The assumption of LTE consists in the following: 
in every "infinitesimal" volume element one can define 
thermodynamic variables in the usual way and these are 
related to each other through the relations which hold 
for the system at equilibrium. It should be emphasized 
from the outset that this assumption is never exact: there 
always exist corrections of the order of the gradients im- 
posed on the system. However, by choosing sufficiently 
small gradients it is always possible to reach a situation 



in which a definition of the temperature and density vari- 
ables as if the volume element was in equilibrium, leads 
to reasonable values for the local chemical potential. Let 
us first discuss the issue of defining the temperature: in 
thermal equilibrium it is well known that the particle 
energies have a Boltzmann distribution. We may there- 
fore define temperature as the parameter characterizing 
this distribution, arguing that temperature is ill-defined 
if the distribution is not Boltzmann. To this end we have 
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FIG. 6: Semi-logarithmic plot of the particle's energy distri- 
bution P X (E) at different positions along a channel obtained 
from the simulation described in Fig. ^ The different curves 
correspond to a fit to a Boltzmann distribution for each posi- 
tion. From these fits, we obtain the temperatures: T = 142.69 
at x/L = 0.1333 (circles), T = 149.70 at x/L = 0.4666 
(squares) and T = 158.69 at x/L = 0.9333 (triangles). The 
curves have been re-scaled for clarity. 
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FIG. 7: Temperature profile of Fig. |B| compared with the tem- 
peratures (open circles) obtained from the fit to the Boltz- 
mann distribution of P(E; x) shown in Fig. H. 



FIG. 9: Window average of the deviation from the Maxwell 
distribution given in terms of the ratio T(v x ) = {P(v x ) — 
P(—v x )} /P(v x ) for a simulation with A(fi/T) = —0.06 and 
T constant in a channel of length L — 30. In the Inset, the 
original result for F(v x ) is shown. 



computed the energy distribution of the gas of particles 
P(E; x) as they cross a narrow slab centered at some po- 
sition x. In Fig. [| we show (in symbols), the results 
for P(E; x) measured at three different positions for the 
same simulation described in Fig. ^| At each position 
the distribution P(E; x) is consistent with the Boltzmann 
distribution, thus indicating that the gas of particles be- 
haves locally as if it was in equilibrium at some tempera- 
ture T{x). If we determine the temperature T(x) by a fit 
of P{E\ x) to the Boltzmann distribution and compare 
these values with the mean energy per particle profile of 
the system shown in Fig. ^, we see in Fig. that T(x) 
coincides with the mean energy of the particles measured 
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FIG. 8: The quantity ln(p/T) (solid line) as a function of 
the position x computed from the density and temperature 
profiles obtained in the simulation described in Fig. |E| The 
dashed line corresponds to a linear profile of the chemical po- 
tential ji/T joining its nominal values in the baths as obtained 
from (h). The solid line has been shifted from its obtained nu- 
merical value for a better comparison. 



locally at the position x. Thus, in the steady state, a lo- 
cal Boltzmann distribution is established. Therefore, the 
identification of the mean energy per particle with the 
local temperature in this model is justified. 

In Fig. |^, we now compare the dependence on x of the 
quantity ln(p/T) with a linear profile between the nom- 
inal values of fi/T in the thermo-chemical baths for the 
same simulation. The agreement between both curves 
(the profile of \n(p/T) has been shifted for comparison), 
indicates that A[i/T = Aln(p/T), supporting the fact 
that the gas of particles inside the channel behaves lo- 
cally as an ideal gas. By this we mean that the relation 
between the chemical potential, the density and the tem- 
perature for an ideal gas in equilibrium holds good locally 
to an excellent approximation in the NSS under study. 

There is, however, an obvious discrepancy: since the 
NSS generally carries a non-zero particle current, it is 
clear that the average velocity is non-vanishing, thereby 
contradicting the Maxwellian distribution for the veloc- 
ities. This is shown in detail in Fig. || where the de- 
viation from the Maxwell distribution given in terms of 
ratio T(v x ) = {P(v x ) - P(-V x )}/P(v x ) is plotted. The 
systematic positive value of T(v x ) is an indication that 
the average velocity is greater than zero |2jJ]. This ve- 
locity, however, is proportional to the particle current, 
and hence to the gradients. Since these must be assumed 
small for LTE to hold, it is a small effect, which vanishes 
in the relevant limit. At this point, it is worthwhile to 
make the following point, when one states that a model 
such as that defined in (TTJ does not satisfy LTE, it means 
that the deviations from LTE do not decrease as the gra- 
dients. In that case, for example, they only decay as the 
imposed temperature difference goes to zero, which in the 
thermodynamic limit is a much more stringent condition. 
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V. NORMAL TRANSPORT AND ONSAGER 
RECIPROCITY RELATIONS 

Having shown that we can assign an unambiguous 
meaning to local thermodynamical quantities in the non 
equilibrium steady states reached by the system, we can 
now study the transport properties of our systems, and 
whether these comply with the predictions of irreversible 
thermodynamics. 

From the general theory of irreversible processes, to 
linear order, the heat and particle currents 3 U and J p 
can be written as follows (see e.g. [p6f): 



L pu V- 



rp L up X7 , 

1 



(12) 



and the Onsager reciprocity relations read in this case 



(13) 
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FIG. 11: Particle's temperature profile for a simulation where 
the quantity fi/T is kept constant. The nominal values for the 
baths' temperature are To = 145 and T\ = 155 in a channel of 
length L — 30. In the inset, the profile of fi/T as calculated 
from (|^) with /Lto = 1 is shown to corroborate its constant 
value along the channel. 



We now wish to show a central feature of our model: 
namely that its transport properties are normal, mea ning 
that the various transport coefficients appearing in fll2| ) 
do not depend on the length of system, and are thus 
well defined in the thermodynamical limit. In Fig. ^ 
we show the dependence of the currents on the length 
L of the system, for a typical realization, in which we 
keep the differences LVT and LVp/T fixed as we vary 
the system size. The 1/L dependence observed confirms 
that transport is normal. 



In order to obtain the value of the coefficients in (12), 
it is enough to perform two simulations: Fixing the value 
of VT and setting V(/Lt/T) =0 yields L uu and L pu from 
the direct measurement of the energy and particle flows, 
while setting VT = and fixing V(/x/T) gives L up and 
Lpp- We have performed simulations with temperature 
and chemical potential differences up to 20% of the min- 
imal nominal values at the walls, and in all cases we have 
found normal transport consistent with ([l3"|). 

In Fig. [ll], we show the particle's temperature and 
chemical potential profiles obtained from a simulation 
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FIG. 10: Size dependence of the heat and matter currents 
for simulations with a fixed temperature difference, AT = 
20, and fj,/T constant. The dotted lines corresponds to 1/L 
scaling. 



with constant /i/T and a temperature difference of AT = 
10 around T = 150 in a channel of length L = 30 aver- 
aged over 470 realizations. After the steady state has 
been reached, we found that both a heat current and a 
particle current were driven by the temperature gradi- 
ent. From ( |l2] ) we obtained the following values for the 
Onsager coefficients: 



L uu = (0.7920 ± 0.0050)pT 5/2 
Lpu = (0.1271 ±0.0017)pT 3/2 



(14) 



From the complementary simulation (see Fig. [12]), 
where temperature is kept constant at T = 150 and a 
chemical potential gradient is imposed with A(/i/T) = 
—0.06 we obtained 



L up = (0.1272 ±0.0048)pT 
L pp = (0.1050 ±0.0030)pT 



.3/2 
1/2 



(15) 



In ( jl4| ) and (|T^) we have written the explicit depen- 
dence of the Onsager coefficients on the particle's den- 
sity p and temperature T. The temperature dependences 
arise from simple dimensional analysis given the fact that 
the system is homogeneous in energy, or equivalently, 
that it has no proper time scale. The same cannot be 
said of the density, but the obtained linear dependence 
shown in Fig. 13L at least in this range of density values, 
is not too surprising. 

The obtained values for the symmetric coefficients L pu 
and L up confirm ([l3]) to within our numerical accuracy. 

As a consistency check, we have also studied a "canon- 
ical" situation, in which we suppressed absorption and 
emission of particles at the walls while still allowing heat 
exchange. In this situation there is no flow of matter in 
the steady state. The relationship between heat flow and 
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FIG. 12: Profile of the quantity /i/T as calculated from (^) 
with fj,o = 1 for a simulation where the nominal values of the 
baths' temperature are set to the same value T = 150. with an 
end-to-end chemical potential difference of A(/j,/T) = —0.06 
in a channel of length L = 30. In the inset, the particle's 
temperature profile is shown. 



temperature gradient becomes J u = kVT, with k given 
by the following expression: 



Luu-L pp L U pL p U 



T 2 L 



(16) 



/>/' 



This relationship was found to hold to good accuracy, 
thus confirming the validity of our "grand canonical" sim- 
ulations by which the L's were evaluated. We remark 
that the coupling between the two currents in the NSS is 
non-trivial in the following sense: in a "canonical" sim- 
ulation, the simplest assumption for the dependence of 
the density on the position is that the trajectory of each 
particle covers the sample uniformly. Then the local tem- 
perature merely determines the speed at which the orbit 
is being traversed. This would imply that in such a situ- 
ation, the particle density should scale inversely with the 
average velocity, that is 



p{x)T{x) 1 ^ 2 = const 



(17) 



In terms of the transport coefficients defined in (12), 
( |l7| ) is equivalent to L pu = (d+l)TL pp /2, where d is the 
dimension of the system. This relation corresponds to a 
system for which all transport arises from uncorrelated 
Markovian motion of the particles, as in the Knudscn 
gas ^0|. However, ( |l7| ) does not hold in our system as 
we always find a small but systematic spatial variation 
in this quantity, the size and sign of which depend on 
the value of rj. Thus, our system cannot be accurately 
described in such a simple manner. 

Despite the fact that the dynamics of the system is 
not ergodic when an external magnetic field is applied, 
we have performed systematic simulations for several val- 
ues of the magnetic field. In particular, we have measured 
both the matter and the energy currents appearing in the 
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FIG. 13: Density dependence of the particle J p and heat J u 
currents. Each symbol in these plots correspond to the av- 
eraged fluxes obtained from a simulation with constant p/T 
and a temperature gradient for which To = 145 and T\ = 155 
in a channel of length L = 30. The lines correspond to linear 
fits. 



direction perpendicular to the applied thermodynamical 
gradients (the so-called Righi-Leduc effect, which is the 
thermal analog of the Hall effect ) • Even when a mag- 
netic field is applied, the existence of the up-down sym- 
metry in the system makes the dynamics equivalent to 
a situation where time-reversal symmetry holds. Due to 
this situation we are not able to verify the validity of the 
Onsager-Casimir relations which are the generalization 
of the Onsager relations when time-reversal symmetry 
is broken. Nevertheless, we have verified that all cross 
L-coefficients satisfy the Onsager relations to within nu- 
merical accuracy. The validity of the Onsager-Casimir 
relations deserves further investigation. 



VI. GREEN-KUBO FORMALISM 

Another interesting question we are able to address 
with our system is whether the Green-Kubo relations 
hold. These relate the phenomenological collective trans- 
port coefficients with the equilibrium time correlation 
functions of microscopic dynamical variables. 

As finite size effects limit the range of validity of the 
Green-Kubo relations, we begin by presenting a deriva- 
tion of these relations which is appropriate to our finite 
length system |^7| . Starting from the phenomenological 
transport equations ( |l2| ) , and using the fact that the gas 
of particles is an ideal gas, we express (|l^) in terms of 
the gradients of the energy density u and particle density 
P- 



- 2L up Vp 



(18) 



T - — ( T Lpu 
p ~ pT \ pp T 



As the analysis is restricted to the linear regime, the 
factors of p and T appearing in the coefficients of the 
gradients are taken as constants. Using energy and mass 
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conservation, we obtain from (US) 



(19) 



where the matrix A is the matrix of coefficients whose 
elements can be readily identified from (|lli|). 

Applying a Fourier transform in space to ( [l9| ) we ob- 
tain 



Mfc(i)\ _ -tfc 2 A f u k (0) 
Pk{t) ' Ufe(O) 



(20) 



where Uk(t) and /3fe(i) are the fc Fourier components of 
u(x, t) and p(x, t) evaluated at time t. Now we take the 
exterior product of ( |20| ) with the vector (u%(0) /5£(0) an d 
average over the ensemble of equilibrium realizations to 
obtain 



C fc (t) 



-tie' A 



u, 



feO 



(21) 



The matrices Cfc(i) and U^o consist of correlation func- 
tions of the densities given by 



C fc (i) 



(fi fc (t)fij(0)) (&*(*)#(<))) 

<M*K(o)> (o)) 



and 



U fc0 = 



_ ( (|fi fc (0)| 2 ) (%(0)pl(0)) 
(p k (0)«J(0)) <|p fe (0)| 2 > 



(22) 



(23) 



The above correlation functions were derived explicitly 
as forward time correlation functions; however, as they 
pertain to the stationary state of the system, we can 
extend them backwards in time as even functions. Having 
done this, the Fourier transform in time of (pl| ) yields 



2fc 2 



A 2 fc 4 



rAUfeo . 



(24) 



If we use again the continuity equations we can express 
the matrix Cfe(u>) in terms of time correlations of the 
energy and mass fluxes as 



fc 2 - 



where the matrix JJ/j is given by 



k[ } ~ VC* JpJ >) c fcJpJ » 



(25) 



(26) 



and C k j a j b {u>) = ^J fea (cj)J fe6 (-a; 

functions in the k-u> space. Inserting ( |25D into (| 
obtain 



are flux correlation 
we 



JfcM 



2lu 2 



A 2 fc 4 



rAUfeo 



(27) 



In the limit k 2 — ► 0, the relation in 
pendent of k and u>: 



A = hk(^ko 



(p7[) becomes inde- 



(28) 



However, for finite fc, the limit attained by taking lo — > 
yields I* (a;) = 0. In ( p8| ) the Onsager coefficients are 
expressed as a function of correlation functions of energy 
and particle fluxes and the static correlation functions of 
the densities. Furthermore, from the equilibrium statis- 
tics of the system for the energy and particle distributions 
of the gas of particles, the matrix Ufeo can be explicitly 
written as 



/27VT 2 NT 
y NT N 



(29) 



where N and T are the number of particles and temper- 
ature of the gas respectively. 

Finally, inserting (|2^) into (E§h we obtain a compact 
expression for the Green-Kubo relations for the Onsager 
coefficients 



L ab - 2jpki a 3 b { W ) 



(30) 



for a, b = {u p}. Here we have used that pL = N. 

In order to obtain a value for the Onsager coefficients 
from the Green-Kubo formulas of (|3^) we measured the 
equilibrium correlation functions of spatial Fourier com- 
ponents of the energy and particle currents. To this ef- 
fect, we performed the measurements in a channel com- 
posed of L = 120 unit cells. These measurements were 
done in microcanonical simulation with periodic bound- 
ary conditions, as it is clear that transport coefficients 
vanish in the static limit if the boundary conditions 
are reflecting. Also, the above relations were obtained 
strictly as the limit of fc — > 0, we therefore worked at 
small finite values for the wave number fc at frequency 
ui, which corresponds to the situation in which (E7j) is 
applicable. Note that the usual way of evaluating the 
diffusion constant by integrating over the current corre- 
lations in the periodic case, which corresponds to having 
lj small and fc strictly equal to zero, does not carry over in 
a straightforward way when the particles are interacting, 
as is the case in our system. 

The periodic microcanonical simulation was performed 
at an energy corresponding to a temperature T = 150 
with N = 658 particles. The particles and discs were 
set up in an equilibrium state, and we allowed the sys- 
tem to further equilibrate by itself before we started the 
measurements. The energy and particle flux correlation 
functions were obtained as follows: We measured at each 
of the n = 1 . . . L cells the energy J nu (t) and particle 
Jnp(i) currents as a function of time. They were mea- 
sured at discrete values of time with At = 0.005 (which 
corresponds to approximately a third of the typical time 
required to traverse the shortest collision path in the sys- 
tem). From these, we obtained the fc Fourier component 
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FIG. 14: Dependence of the different flux-correlation func- 
tions on ui. The solid lines correspond to Cfej a j 6 (a;) for k = 1. 
The dashed lines are for k = 10. The dependence on omega 
of Cfej u j p (ct;) is the same that for Cfcj p j u (a;) within the nu- 
merical accuracy and therefore, not shown. 



of the fluxes as 

L 

J kp (t) = J2 el27rkn/L ^ P (t) (31a) 

n=l 
L 

hu(t) = X! ei27rfe " /ij ««(*) • ( 31b ) 



From a time series of 4x 10 7 data, approximately equiv- 
alent to 10 6 times the mean free flight time, we obtained 
the time averaged Fourier component of the energy and 
particle flux time-correlation functions 



CfcJ.J^f) = <J*a(*)J^(0)) t 



(32) 



Finally, we performed the time Fourier transform of these 
quantities. In Fig. |lj we show the plots of the correlation 
functions Cfcj a j b (w), as functions of u> for values of k = 1 
and 10. These plots show that at low frequencies the cor- 
relations functions tend to zero, in agreement with (27). 



However, there is a clear plateau in the plots correspond- 
ing to the regime in which the Green-Kubo relations hold. 
At higher frequencies, the curves again deviate from the 
plateau. This is due to the fact that these high frequen- 
cies correspond to times that are shorter than typical 
collision times, for which the diffusive transport assump- 
tions used to derive (|27]) no longer hold. Taking the 
plateau values for the correlation functions, which are 
essentially independent of k and to, we can compute the 
Green-Kubo predictions for the Onsager coefficients us- 
ing (p0|). In Table || we summarize the values obtained 
for the different Onsager coefficients and compare them 
with those obtained previously at this temperature and 
density. 

An interesting issue involving the linear-response of the 
system is the macroscopic equivalence between transport 





Green-Kubo 


Gradients 


Lpp 


0.1050 ±0.003 


0.1030 ±0.002 


L pu 


0.1276 ±0.0016 


0.1271 ±0.0017 


Lup 


0.1276 ±0.0016 


0.1272 ± 0.0048 


Luu 


0.7920 ± 0.012 


0.7710 ±0.005 



TABLE I: Onsager coefficients calculated from the Green- 
Kubo formulas ( pOj ) compared with those values obtained pre- 
viously from direct measurements. 



processes driven by thermodynamical forces (due to gra- 
dients in thermodynamical quantities) and those driven 
by real forces (like externally applied fields). 

It has been argued by van Kampen |2^| that the lin- 
ear response theory derivation of Green-Kubo formula for 
the electrical conductivity is not correct as the response 
of the microscopic trajectories to an applied electric field 
can not be taken as linear (see also J2j| ||l| ) . It has been 
argued by Visscher ]2j|, however, that a difference may 
exist between the case of mechanical forces, for which 
van Kampen's objection might hold, and thermodynam- 
ical forces such as temperature and chemical potential 
gradients, for which linear response should provide cor- 
rect answers. In order to test this point, we performed 
simulations for our system in an electric field (all par- 
ticles having the same charge and no interaction) and 
compared this with the result of the corresponding gra- 
dient in chemical potential. 

At the level of the Green-Kubo relations the equiva- 
lence follows directly from the linear response results for 
a chemical potential gradient or an applied electric field: 
Both are related to the autocorrelation function of the 
particle current when all particles have the same charge, 
which is the case in our system. 

In order to test these ideas we have performed two com- 
plementary simulations: In a first simulation we impose 
a chemical potential gradient at constant temperature. 
In a second simulation both temperature and chemical 
potential are constant and an external uniform electric 
field is applied to the system in which the particles now 
carry a unit electric charge q. To compare the macro- 
scopic flows obtained in these simulations, the magnitude 
of the electric field E is fixed to 



A^i = qEL , 



(33) 



so that, the work done by the electric field in taking a 
particle from one side of the channel to the other is the 
same as the chemical potential difference. In (|33|), L is 
the size of the system. 

In Fig. the temperature profile of both simulations 
is compared. When the electric field is applied (Fig. [l5|- 
a) , the temperature in the bulk of the system increases 
due to the internal dissipation of the work done by the 
field. This is the well known Joule heating effect [ p6[ . 
In this case, we have checked that the dependence of the 
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FIG. 15: Temperature profiles obtained from a simulation at 
constant T — 150 and mean number of particles N ~ 41.2 
in a channel of length L = 30 obtained a) from a simulation 
with constant fi and a applied electric field of magnitude E = 
0.5. and b) without electric field and an imposed chemical 
potential gradient of V(/Li/T) = 0.1. The circles correspond 
to the averaged kinetic energy of the discs for the respective 
simulation. The dashed line is the expected value for the 
profile. 



increment of temperature in the bulk AT is quadratic 
in the field and observed that it holds even far beyond 
the linear regime. In the case of an imposed chemical 
potential gradient (Fig. |l^-b), the temperature profile 
appears to be flat, indicating that in this case no effect 
equivalent to Joule heating is present. 

We have also measured the particle density profile and 
obtained the profile for the quantity [ijT computed us ing 
the relation (Q) in both cases. These are shown in Fi g, plf - 
a for the case of an applied electric field and in Fig. |16|-b 
for an imposed chemical potential gradient. In case(a) 
the chemical potential is flat since we are not taking into 
account the electric contribution. In case (b) the profile 
is linear and shows quadratic deviations from linearity 
similar to those observed for the temperature profile in 
case (a). This is shown in the inset of Fig. |l6|-b where 
we subtracted the theoretical linear /i/T profile. 

Finally, for the heat and matter flows we have obtained 
for the case of an applied electric field 



-0.00829 ±0.00006 
-1.59 ±0.02 , 



(34) 



and for the case of an imposed chemical potential gradi- 
ent 



J p = -0.00829 ±0.00008 
J„, = -1.59 ±0.04 . 



(35) 



Thus, the fluxes obtained in both simulations corroborate 
our initial expectations within the numerical accuracy. 



VII. RATES OF ENTROPY PRODUCTION 
AND PHASE SPACE VOLUME CONTRACTION 



FIG. 16: Profiles for the quantity fi/T minus the expected 
profile /io/T obtained from the simulations described in 
Fig. ^ with (a) an applied electric field and (b) with an im- 
posed chemical potential gradient. 



eral, this rate has been at the focus of considerable in- 
terest. In particular, in recent work on thermostated 
systems, the average rate of phase space volume contrac- 
tion has been identified with the rate of entropy produc- 
tion [§ §, @, §, ||. While this work is certainly 
of considerable interest, it is vital to be able to under- 
stand how such claims might generalize to purely Hamil- 
tonian systems, for which phase space volume is rigor- 
ously conserved, at least in the fine-grained sense, due 
to Liouville's theorem. Our system is well-suited to this 
purpose, since its dynamics is quite transparent and its 
equilibrium behaviour is trivial. 

In our system, it is clear that phase space volume can 
only be created or annihilated at the boundaries, due to 
the conservation of phase space throughout the internal 
dynamics (note that the collision rules (||) are volume 
preserving). To estimate the variations in phase space 
volume induced by the incoming and outgoing particles, 
we proceed as follows: let us denote by T\ and T 2 the 
two temperatures and by /ii and /12 the two chemical po- 
tentials at either end of the system (once again, due to 
LTE there is a local "temperature" and "chemical poten- 
tial" in the stationary state of the system). A particle 
entering at one end of the system therefore contributes 
to the change in phase space volume f2 by changing both 
the the local energy per unit volume and the local par- 
ticle density. Since, as follows from the assumption of 
LTE, the local thermodynamical variables are connected 
to each other in the usual way, the increase in phase space 
volume is given by 



Alnn = g(i(AE0 1 -|(AJV) 1 ) 



(36) 



where i can take the values 1 and 2 for either end. Due 
to stationarity, however, it is clear that 



It is generally argued that stationary states are char- 
acterized by minimal rate of entropy production. In gen- 



(AN), 
(AU\ 



(AA0 2 
(At/) 2 



= j«At 



(37) 
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Putting ( p7| ) into (p6|), one obtains the expression for 
the rate of production of the phase space volume: 



dt 



L 



j«V-+j p V 



(38) 



Using (|12|), and the positive definiteness of the matrix 
of transport coefficients L a p, one finds that the r.h.s of 
( |38| ) is always negative: on the whole, phase space vol- 
ume flows out of the system into the baths. This can 
be interpreted in two complementary ways. At first, one 
might argue that if a system has ever contracting phase 
space volume, it must eventually end up on an invariant 
set of measure zero. This is in fact quite reasonable and 
in good agreement with various findings on the dynamics 
of reversible thermostated systems. Indeed, various peo- 
ple JyJ |32], |33|, Q have found that the invariant mea- 
sure for such systems is singular with respect to Liouville 
measure. If we were to generalize such an assumption to 
our system, it would mean that the invariant measure of 
the stationary state it eventually reaches is also singular 
with respect to Liouville measure, which would in turn 
signify that the total phase space volume of the steady 
state is zero. In principle, this seems amenable to a nu- 
merical test, but this is not really the case: due to the 
rather high dimensionality of our system, the fractality 
of the support of the steady state measure is presum- 
ably extremely difficult to observe, if it is at all possible. 
Therefore, even if we could connect entropy production, 
say, with some kind of escape rate (as proposed by Gas- 
pard and coworkers |p5|), it is not at all clear how to 
measure such a quantity in even so simple a model as 
ours. Whether truly low-dimensional systems exist, for 
which such quantities are readily accessible and which 
also exhibit LTE and normal transport, is an interesting 
open problem. 

Another interpretation of (|3^) is the following: in the 
stationary state, because we are in fact in LTE, that is, 
approximately in a state of thermal equilibrium in every 
volume element, we may think that at least the coarse 
grained phase space volume, measured by the thermo- 
dynamical entropy, will remain constant in steady state. 
But if, as claimed in (|38|), an amount of phase space vol- 
ume leaves the system each unit of time, then an equiva- 
lent amount of coarse-grained phase space volume must 
be produced inside the system. That this last statement 
does not contradict Liouville's theorem is, of course, well- 
known from traditional examples in statistical mechanics. 
Such an interpretation might then justify the interpre- 
tation of the r.h.s of (38) as minus the rate of entropy 
production, in conformity with the results stated for re- 
versible thermostated systems. 

To make such an identification plausible, we need a def- 
inition of the entropy out of equilibrium. It is well-known 
that in the general case such a definition poses formidable 
problems. However, since our system is an ideal gas but 
little perturbed from equilibrium, it is possible to use 
Boltzmann's expression for an entropy density per unit 
volume in terms of the one-particle distribution function 



f(x,v;t), defined by 



s(x) 



dv f(x,v;t)lnf(x,v;t) 



(39) 



From this one may give an elementary expression for the 
rate of entropy production through the following consid- 
erations. Due to LTE, the stationary distribution /(a;, v) 
is given by a local Maxwellian 

f(x, v) = p{x)f}{x) x ' 2 exp [-l3(x)v 2 /2] + 5f(x, v) (40) 

where Sf(x, v) is a correction term of the order of the im- 
posed gradients, which accounts for the currents in the 
system. As shown in (|9|) such a correction has been ob- 
served in our system. It then seems reasonable to divide 
the time variation of f(x, v) into a convective part given 
by —vf x plus a collisional part, of which we need say 
nothing. In this case, we may define the entropy current 
density as follows 



.Is = - / dvvf(x,v)lnf(x,v) (41) 

= -J dvvSf(x,v)[\nf(x,v) + 1] . (42) 

(43) 

Inserting @) into (ph yields 



M. 1. 



(44) 



In the stationary state, the local rate of entropy pro- 
duction is given by the divergence of the entropy current. 
If one takes it as given by ( f44| ) , we obtain again the rate 
of entropy production given by the standard expression, 
this time with the correct positive sign. 



VIII. CONCLUSION 

We have introduced a simple model system character- 
ized by reversible Hamiltonian microscopic dynamics, the 
equilibrium properties of which are those of an ideal gas 
when viewed as a thermodynamic system. We performed 
extensive numerical studies of the system both in equilib- 
rium and in non-equilibrium steady state. We find that 
when driven out of equilibrium, the properties of the sys- 
tem are consistent with the hypothesis of local thermal 
equilibrium, that is, in the stationary state, one finds 
a local Boltzmann distribution for the energies, leading 
therefore to an unambiguous definition of the local tem- 
perature, chemical potential, etc. In this situation its 
transport properties are found to be entirely similar to 
those of realistic interacting many particle systems: it ex- 
hibits coupled mass and heat transport and the two cross 
transport coefficients satisfy Onsager's relations. When 
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time reversal invariance is broken by means of an ap- 
plied constant magnetic field, the appropriate generaliza- 
tions of the Onsager relations hold. This is the case even 
though this system turns out not to be ergodic. Further- 
more, we find extremely satisfactory agreement between 
the values of the transport coefficients as deduced from 
the equilibrium dynamics via the Green-Kubo formulas 
and those obtained via direct simulation. A final result of 
interest is that we were able to show equivalence between 
an applied electric field and a chemical potential gradi- 
ent. This result, while a straightforward consequence of 
linear response, does not appear at all obvious from the 
microscopic point of view. It may therefore be considered 
as a further confirmation of the validity of the linear re- 
sponse formalism 

The fact that all these features are present in such 
a simple model, suggests that it may serve as an ideal 
framework to gain insight of how macroscopic transport 
phenomena arise in real systems. The model is also 
well suited to test theories for the description of systems 
in out of equilibrium states. Among these, a study of 
the Cohen-Gallavotti theorem might be very interesting. 
Also, issues linked to those involving entropy production 
discussed in the section VII, in particular the question 



the master equation 



of characterizing the invariant measure of the stationary 
state appears very challenging. 

Finally, extensions to more complicated systems can 
also be considered. By varying the value of the moment 
of inertia of the rotors one can study transport in het- 
erogeneous structures, such as junctions, layered systems, 
etc. Other extensions could include mixtures of scattered 
particles with different masses. For all these extensions, 
most if not all of the equilibrium properties are straight 
forward, and theories for the phenomenological behaviour 
of such systems can be carefully tested. 
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APPENDIX A: VARIOUS POSSIBLE 
THERMO-CHEMICAL BATHS 



In this appendix we present a short discussion of 
Markov processes that can be used to generate a canoni- 
cal or grand canonical ensemble. A fairly general Markov 
process combined with a Hamiltonian dynamics satisfies 



dtPN(x) 



{H N ,p N } + J2 ! W{N'N-x'x) PN ,{x l )dx l 

N' ^ 



~Pn{x) 



W{NN';xx')dx'. 



(Al) 



Here N is the number of particles, x is an abbreviated 
notation for the vector (p, q) and pn{x) is the probability 
density of being at the phase point x having altogether 
N particles. The W(NN'; xx') are the transition rates of 
going from a phase point x with N particles to a phase 
point x' with N' particles. If we are looking for rates such 
that they generate, say, the grand canonical ensemble 
with a given temperature and chemical potential, it is 
sufficient that they satisfy the detailed balance condition 



W(NN';xx') 
W{N'N-x'x) 



= exp[Pfi(N - N') ~(3{H N> (x r ) ~ H N (x))\ . 

(A2) 

Let us now assume that the stochastic part of the dynam- 
ics is limited to the case in which x belongs to a small 
subset r of the full phase space, otherwise the dynam- 
ics is purely Hamiltonian. This corresponds to the case 
of stochastic walls treated in this paper. Assume further 
that, as in our model, the Hamiltonian in T is one-particle 
only, say a pure kinetic energy term. Assume finally that 
the only changes we consider will be the introduction and 
destruction of a single particle. Then the following rates 
are a solution of QA2|): 

W r (N,N + 1) = Xe-^+ri Wr(N + l,N) = A 

(A3) 

If we want this stochastic process to act as a localized 
bath in the system, we take the region T in phase space 
defined by the condition that at least one particle is in the 
region V of the (one-particle) configuration space. The 
rates in (A3) then mean that any particle entering V is 
annihilated at a rate A. Further, particles of speed v are 

being created inside V at a rate \e~^^ + ^ . Since such 
particles may be immediately reabsorbed, however, the 
problem of the particle flux emitted by V is not quite 
straightforward. If V is narrow (of thickness I) in one 
dimension around a given hypersurface S, however, and A 
is ultimately made to go to infinity as I — > 0, the problem 
can be solved as follows. For simplicity assume that all 
momenta always point in the direcion of one of the sides 
of V. The probability that the particle will come out at 
all, conditioned on its having been created at a distance 
r of the side from which it must leave V is given by 



p(r) = cxp [-Xr/v n ] , 



(A4) 



where v n is the velocity normal to the surface S. The 
flux of particles coming out of V is then given by 



Ae"^ 2 / 2 j drp(r) = v n e~^' 2+ ^ (l - e^") , 

(A5) 
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where a is the limiting value of XI. In the limit in which 
ol — ► oo, this model reduces to the one described in the 
text, which was the one primarily used in the simulations 
presented in this paper. The above model can be gen- 
eralized in different ways: in particular, we may allow 
for thermalization without change in the number of par- 
ticles. A possible candidate for such a reaction rate is 
given by 

W(NN; xx') = X' exp [-(3{v 12 - v 2 )/2] 6(q - q'), (A6) 

where again the process acts only in the region T de- 
scribed above. This corresponds to a process in which 
the particle velocities are thermalized while leaving their 
positions fixed, which is, as is readily seen, the algorithm 
for thermalization through collisions introduced in the 
text. In particul ar, t he case a = and a finite rate of 
the type given is (A6) leads to the algorithm we used for 



generating the canonical ensemble. 



Various simulations were made with finite values of a 
in order to test the method. In equilibrium it was al- 
ways found that the correct nominal temperatures and 
chemical potentials were attained after sufficient time. 
In non-equilibrium situations, however, for small values 
of a, the thermodynamic parameters of the system can 
differ from the values imposed by the baths. In partic- 
ular, the energy density showed jumps localized in the 
vicinity of the walls. These energy gaps have been fre- 
quently observed in simulations of other transport mod- 
els 0, H ||, and were studied in §9|. In our 
model, they eventually disappear as we increase the value 
of a, meaning that the absorption of the wall increases 
and the particles (which thermalize the system with the 
bath) are exchanged more easily. 
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